clc;
syms x y;
syms a positive;
z=a-x-y;
f=x*y*z;


I=surf_integral(f,z,[x,y],[0,a-x],[0,a])